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ABSTRACT 

We have simulated the formation and evolution of large-scale structure in the universe, for 
68 different COi3£^normalized cosmological models. For each cosmological model, we have per- 
formed between 1 and 3 simulations, for a total of 160 simulations. This constitutes the largest 
database of cosmological simulations ever assembled, and the largest cosmological parameter 
space ever covered by such simulations. We are making this database available to the astronom- 
ical community. We provide instructions for accessing the database and for converting the data 
from computational units to physical units. 

The database includes Tilted Cold Dark Matter (TCDM) models. Tilted Open Cold Dark 
Matter (TOCDM) models, and Tilted A Cold Dark Matter (TACDM) models. (For several sim- 
ulations, the primordial exponent n of the power spectrum is near unity, hence these simulations 
can be considered as "untilted.") The simulations cover a 4-dimensional cosmological parameter 
phase space, the parameters being the present density parameter fio, cosmological constant Aq, 
and Hubble constant Hq, and the rms density fluctuation as at scale Sh^^Mpc. All simulations 
were performed using a P'^M algorithm with 64^ particles on a 128'^ mesh, in a cubic volume 
of comoving size 128 Mpc. Each simulation starts at a redshift of 24, and is carried up to the 
present. More simulations will be added to the database in the future. 

We have performed a limited amount of data reduction and analysis of the final states of 
the simulations. We computed the rms density fluctuation, the 2-point correlation function, the 
velocity moments, and the properties of clusters. The details of these calculations are presented 
in the full version of this paper. In this abbreviated version, we present only the analysis of the 
rms density fluctuation and two-point correlation function. Our results are the following: 

(1) The numerical value cTg^™ of the rms density fluctuation differs from the value cr™"* 
obtained by integrating the power spectrum at early times and extrapolating linearly up the the 
present. This results from the combined effects of discreteness in the numerical representation 
of the power spectrum, the presence of a Gaussian factor in the initial conditions, and late-time 
nonlinear evolution. The first of these three effects is negligible. The second and third are 
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comparable, and can both modify the value of a% by up to 10%. Nonlinear effects, however, are 
important only for models with erg > 0.6, and can result in either an increase or a decrease in erg. 

(2) The observed galaxy two-point correlation function is well reproduced by models with 
cTg ^ 0.8, nearly independently of the values of the other parameters, fio, Aq, and Hq. For models 
with cTg > 0.8, the correlation function is too large and its slope is too steep. For models with 
CTg < 0.8, the correlation function is too small, its slope is too shallow, and it often has a kink at 
separations of order 1 — 3 Mpc. 

(3) At small separations, r < 1 Mpc, the velocity moments indicate that small clusters have 
reached virial equilibrium, while still accreting matter from the field. The velocity moments 
depend essentially upon Q.q and cxg, and not Aq and Hq. The pairwise particle velocity dispersions 
are much larger than the observed pairwise galaxy velocity dispersion, for nearly all models. 
Velocity bias between galaxies and dark matter is needed to reconciled the simulations with 
observations. 

(4) The cluster multiplicity function is decreasing for models with cjg ~ 0.3. It has a horizontal 
plateau for models with erg in the range 0.4—0.9. For models with ag > 0.9, it has a U shape, which 
is probably a numerical artifact caused by the finite number of particles used in the simulations. 
For all models, clusters have densities in the range 100-1000 times the mean background density, 
the spin parameters A are in the range 0.008 — 0.2, with the median near 0.05, and about 2/3 of 
the clusters are prolate. Rotationally supported disks do not form in these simulations. 

Subject headings: cosmology: theory — large-scale structure of the universe - methods: numerical 
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1. INTRODUCTION 

1.1. Importance of Numerical Simulations in Cosmology 

Observations of the nearby universe reveal the existence of the large-scale structure. The visible matter 
is clumped into galaxies, and these galaxies are not distributed uniformly into space, but instead grouped 
into structures such as clusters, filaments, and walls, separated by deep voids. Velocity structures (deviations 
from Hubble flow) arc observed as well. By contrast, observations of the Cosmic Microwave BackgroLmd 
(CMB) reveal that the universe was extremely uniform near the epoch of recombination. Hence, the present 
large-scale structure must result from an evolutionary process that took place between recombination and 
the present. The most widely accepted scenario assumes that the present large-scale structure originates 
from the growth, by gravitational instability, of primordial density fluctuations present in the early universe. 
Any fluctuation larger than the Jeans length can grow by gravitational instability once the universe becomes 
matter-dominated. If the primordial density fluctuations originates from a Gaussian random process (the 
usual assumption), then the primordial density field is entirely described in terms of its power spectrum 
P{k). The particular form of the power spectrum essentially depends upon the amount and nature of the 
various components (baryonic matter, dark matter, cosmological constant, and so on) present in the universe. 
If we assume a certain power spectrum, we can describe the primordial density field, and the formation 
and evolution of large scale-structure in the universe becomes an initial value problem: starting from the 
primordial density field, we can compute its evolution using the laws of general relativity. Unfortunately, 
this initial value problem is far too complex to be solved analytically. We can simplify the problem by 
noticing that the largest structures observed in the universe are significantly smaller than the horizon. This 
enables us to describe the evolution of the large-scale structure using Newtonian mechanics instead of general 
relativity (Peebles 1980, Chapter 2). Even so, the general problem cannot be solved analytically. This leaves 
two possible approaches: analytical approximations, or numerical simulations. 

Two different kinds of analytical approximations have been considered. The first one is based on the 
fact that the initial fluctuations are small. We can expand the equations describing the evolution of these 
fluctuations in powers of the fluctuations, and solve them using perturbation theory. This approach is 
extremely useful in describing the early evolution of the fluctuations, and has led to very important results. 
However, it becomes inapplicable as soon as the fluctuations become nonlinear. Such fluctuations still have 
to grow by a factor of 10^ to reach the density of a cluster of galaxies, and 10^ to reach the density of a galaxy. 
Clearly, perturbation theory cannot be used to describe the late stages of large-scale structure formation. 
The second analytical approach consists of considering systems with a particular geometry (see, for instance, 
Zel'dovich 1970; Peebles 1980, §§19-21; Fillmore & Goldreich 1984a, 1984b; Bertschinger 1985a, 1985b). The 
most popular analytical models for large-scale structure formation are the Spherical Model, which assiimcs 
spherical symmetry, and the Pancake Model, which assumes planar symmetry. An important assumption of 
some of these models is that the system considered is isolated. For instance, the spherical model can describe 
the evolution of a self-gravitating spherical overdensity, but we must assume that any tidal influence from 
nearby structures can be neglected, an assumption that might be valid at late times but certainly not at 
early times. 

These analytical approximations can therefore describe the universe at early times or at late times, but 
not both. This problem can be solved by using mixed schemes, that combine various analytical approxima- 
tions in a way that allows an analytical description of the evolution of large-scale structure at all epochs. 

The most important ones arc the Press-Schechter Approximation (Press & Schechter 1979), which combines 
perturbation theory with the spherical model, and the Zel'dovich Approximation (Zel'dovich 1970), which 
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combines perturbation theory with the pancake model. 

The alternative consists of using numerical simulation. Unlike analytical models, numerical simulations 
suffer from problems such as limited resolution and numerical noise. Also, simulations provide very little 
insight into the physical processes taking place, compared to analytical models. However, numerical simula- 
tions can describe the evolution of the large-scale structure entirely, from the initial conditions all the way to 
the present, without making any approximation or imposing any restriction on the geometry of the system. 
Cosmological N-body simulations have played a central role in the study of the formation and evolution of 
large-scale structure in the universe during most of the last two decades. These simulations have contributed 
the improve our understanding of the physical process of gravitational instability that leads to structure 
formation, have enable us to conceive and test various cosmological scenarios, and have produced simulated 
universes that can directly be compared with observations (Efstathiou & Eastwood 1981; Centrella & Melott 
1983; Klypin & Shandarin 1983; Miller 1983; Shapiro, Struck-Marcell, & Melott 1983; White, Frenk, & Davis 
1983; Davis et al. 1985; Efstathiou et al. 1985; Barnes & Hut 1986, 1989; Evrard 1986, 1987; Melott 1986; 
White et al. 1987a, 1987b; Frenk et al. 1988; Gramann 1988; Carlberg & Couchman 1989; Villumsen 1989; 
West, Oemler, & Dekel 1989; Couchman 1991; Fukushige et al. 1991; Hernquist, Bouchet, & Suto 1991; Mar- 
tel 1991a; Moutarde et al. 1991; West, Villumsen, & Dekel 1991; Bouchet & Hernquist 1992; Fry, Melott, & 
Shandarin 1992; Park et al. 1992; Bahcall, Cen, & Gramann, 1993; Gramann, Gen, & BahcaU 1993; Melott 
& Shandarin 1993; Babul et al 1994; Pen 1995; Golombi, Bouchet, & Hernquist 1996; Moore, Katz, & Lake 
1996: Yess & Shandarin 1996; Klypin, Nolthenius, & Primack 1997; Kravtsov, Klypin, & Khokhlov 1997; 
Navarro, Frenk, & White 1997; Gross et al. 1998; Thomas et al. 1998) .0 



1.2. The Standard Model 

Theoretical developments in particle theory and early-universe physics, combined with numerical simu- 
lations and observations of the large-scale structure of the universe, led to the emergence during the 1980's of 
what became known as the Standard Cosmological Model. The inflationary scenario (Guth 1981; Linde 1982; 
Albretch & Steinhardt 1982) requires that the universe is spatially flat to a great accuracy. In a matter- 
dominated universe, in the absence of any exotic components such as a nonzero cosmological constant A, 
this requires that the mean density of the universe is equal to its critical density, or, alternatively, 

flo^l, (1) 

where = SttGpq/SHq is the density parameter, po is the mean density of the universe, and Hq is the 
Hubble constant (throughout this paper, we use subscripts to designate present values). In this scenario, 
the large-scale structure of the universe that we observe today results from the growth, by gravitational 
instability, of small density perturbations present at recombination, which originate from quantum processes 
in the early universe. There are two difficulties with this scenario. First, there is strong evidence that the 
amount of "ordinary matter" in the universe is insufficient to satisfy equation (1). Primordial nucleosynthesis 
provides stringent upper limit to the baryonic content of the universe, and shows that the present baryon 
contribution to the density parameter, JIboj is less than 0.026/i~^, where h is the Hubble constant in units 



''Cosmological numerical simulations have been used intensively since 1981, and their results have appeared in hundreds of 
publications, so this list is necessarily incomplete. We decided to include only the key publications by each research group. We 
also excluded one- and two-dimensional simulations (for brevity), and simulations with hydrodynamics, which involve the next 
generation of numerical algorithms. 
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of 100kms~^Mpc~^ (Krauss & Kernan 1995; Copi, Schramm, & Turner 1995; Krauss 1998). Furthermore, 
dynamical studies of rich chisters of galaxies show that the contribution to the density parameter of the 
matter that clusters at that scale is ficiusters 

= 0.2 ± 0.1 (Gott et al. 1974; Carlberg et al. 1996; Lin 
et al. 1996). Second, observations of the temperature fluctuations in the Cosmic Microwave Background 

(CMB) provide an upper limit to the density fluctuations at recombination, of order Spj p ^ 10~^. Such 
fluctuations would grow by gravitational instability, to reach an amplitude of order 10~^ by the present, 
clearly insuflicient to explain the origin of large-scale structure and galaxies. 

These two difficulties were solved by postulating the existence of a component known as dark matter. 
This dark matter, which can be detected only through its gravitational influence, makes up the difference 
between the amount of matter required to satisfy equation (1) and the amount of matter that is observed 
or indirectly measured. We can then reconcile the dynamical estimates of SIq with equation (1) by assuming 
that the dark matter is distributed more smoothly than the luminous matter that makes up galaxies and 
clusters, an idea known as biasing (Kaiser 1984). However, since the dynamical estimates of VLq exceed the 
limit imposed by primordial nucleosynthesis, some amount of dark matter must be clustered on galactic and 
cluster scales, even though the bulk of the dark matter is smoothly distributed into space. Also, density 
fluctuations in the dark matter start growing when the universe becomes matter-dominated, and by the 
time recombination occurs, the dark matter fluctuations have already grown by a factor of ^ 2Qh? 
(Kolb & Turner 1990, §9.5). These fluctuations provide potential wells into which the baryons fall soon 
after recombination. This enables us to reconcile the CMB measurements, — which are sensitive to the 
fluctuations in the baryonic matter, but not the dark matter — with the existence of galaxies and large-scale 
structure. Many candidates for dark matter have been suggested, several of them emerging from particle 
theory. These various forms of dark matter are usually classifled as Hot Dark Matter (HDM) and Cold Dark 
Matter (CDM). 

Cosmological numerical simulations performed in the 1980's have shown that an = 1 universe in which 
most matter is in the form of Cold Dark Matter, and the distribution of luminous matter is biased relative 
to the distribution of dark matter, can successfully reproduce all observations of large-scale structure that 
were then known (see, e.g. Davis et al. 1985), while satisfying the constraints imposed by the inflationary 
scenario, primordial nucleosynthesis, and the CMB (unlike HDM, which as difficulties explaining structure 
formation at galactic scales). These simulations played a central role in establishing the Oq = 1, biased CDM 
model as the Standard Cosmological Model. 

1.3. Non-standard Models 

The Standard Model, which has been hailed by many as the "final answer" to the problem of structure 
formation and evolution in the universe, ran into serious problems during the 1990's. In this subsection, we 
brieffy review these problems. 

1.3.1. The Age Problem 

In a flat, matter-dominated universe, the age of the universe to is 2/3 of the Hubble time, that is 
to = 2/3ifo = 6.52 X 10^/i~^years. For h in the range 0.5 — 1, this corresponds to an age in the range 
6.52 — 13.04 X 10^ years. Measurements of the ages of globular clusters indicate that the oldest clusters are 
certainly older than 9.5 x 10^ years, and most likely in the range 11 — 13 x 10^ years (Jimenez et al. 1996; 
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Chaboyer 1998; Chaboyer et al. 1998; Jimenez 1998). These measurements are only marginally consistent 
with the standard model, in the sense that they require a Hubble constant near its smallest possible value, 
h ~ 0.5. However, recent observations have significantly reduced the range of plausible values for the Hubble 
constant, showing that h is likely to be in the range 0.65 — 0.75 (Freedman, W. L. 1998, and references 
therein). In the standard model, this corresponds to an age in the range 8.69 — 10.03 x 10^ years. The upper 
end of this range is still compatible with the measured ages of globular clusters, but just barely. 

1.3.2. Large- Scale Structure 

Until 1992, the amplitude of the primordial density fluctuation power spectrum was unknown. We were 
free to tune this amplitude in order to reproduce the correct amount of galaxy clustering observed today. 
The latter is usually characterized by the rms density fluctuation erg at a scale of 8/i^^Mpc. Observations 
of clusters of galaxies show that erg ~ 0.6 for the standard model (Viana & Liddle 1996, see eq. [42] below). 
The discovery by the COBE DMR experiment of degree-scale fluctuations in the CMB temperature (Smoot 
ct al. 1992) has eliminated this freedom, by fixing the amplitude power spectrum. A COi3i?-normalizcd 
Standard Model produces too much structure at cluster scales (Barlett & Silk 1993). The resulting value of 
(Tg is 1.22 for h = 0.5 (Bunn & White 1997), too large by a factor of 2, and becomes even larger for larger h. 
By combining the COBE result with observations of the present large-scale structure, we obtain a constraint 
on the quantity Q.oh, which is 

0.2 < ^loh < 0.3 (2) 
(Peacock & Dodds 1994). Since h is certainly larger than 0.5, this implies fio < 0.6. 

1.3.3. The Baryon Catastrophe 

In the Standard Model, the baryon fraction of the universe is small. Primordial nucleosynthesis imposes 
the constraint that Obq/i^ < 0.026 (Krauss & Kernan 1995; Copi, Schramm, & Turner 1995; Krauss 1998). 
For h = 0.65, this corresponds to flso = 0.061. Hence, if = 1, at most 6% of the matter if the universe 
is composed of baryons, the rest being dark matter. However, observations of X-ray clusters reveal that the 
baryon fraction in these clusters is ~ O.l/i"^'^, or 19% for h = 0.65 (Briel, Henry, & Boringer 1992). Hence, 
X-ray clusters contain a large excess of baryons relative to dark matter compared with the average values in 
the universe, a situation often referred to as the baryon catastrophe. This problem could be solved if we can 
think of a physical process that would concentrate the baryons inside clusters, creating a bias relative to the 
dark matter. However, no such physical process is known. It is much simplier to assume that the density 
parameter is less than unity. In this case, the universal baryon fraction is not CIbo, but Qbo/^o- This baryon 
fraction is ^ 0.1h~^'^ according to observations of X-ray clusters, and smaller than 0.026/ClQh^ according to 
primordial nucleosynthesis. Combining these two results, we get 

Qoh^^'^ < 0.26 . (3) 
This rules out the Standard Model (unless h < 0.07 !). For h > 0.5, this limit becomes f^o < 0.37. 
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1.3.4- Evolution of Cluster Abundance 

In the Standard Model, the density parameter fl is unity at all times, and density fluctuations can grow 
by gravitational instability all the way to the present. In other models, density fluctuations can grow at early 
times, when is near unity. But eventually, Q drops significantly below unity, and the density fluctuations 
"freezes-out." Hence, in a model with ^ig < I, the present abundance of clusters should be comparable to 
the abundance immediately after freeze-out, since not much growth has taken place since then. Conversely, 
in the Standard Model, the present abundance of clusters should be larger than the past one, since the 
growth of density fluctuations never freezes out. Bahcall, Fan, & Cen (1997) and Bahcall & Fan (1998) have 
determined the mass of three massive distant clusters, located at redshift z > 0.5, and showed that in a Oq 
universe, there should be only ~ 10~^ such clusters at 2: > 0.5. They conclude that the density parameter 
is in the range 

0.1 < 1^0 < 0.35 (4) 

(BahcaU 1999). 

1.3.5. Distant Type I Supemovae 

The relationship between the Imninosity distance Dl and the redshift is model-dependent. If standard 
candles can be observed at cosmological distances, then the Dl{z) relationship can be inferred, and limits 
can be placed on the value of the cosmological parameters. This method was recently applied to samples 
of distant ("High-z") Type I supcrnovae, by two independent research teams (Garnevich et al. 1998, and 
references therein; Perlmutter et al. 1998, and references therein). Applied to models with a nonzero 
cosmological constant A, their observations provide severe constraints in the Oq — ^0 phase space (where 
Ao = A/'SH^). Not only the Standard Model is excluded with a high degree of confidence, but open models 
(Oq < 1)7 without a cosmological constant are also excluded, unless ^Iq is very small. Observations of the 
CMB (White 1998; Tegmark et al. 1998) provide a different constraint, that does not rule out the Standard 
Model. However, combining the CMB and Type I supernovae observations leads to separate determinations 
of Oq and Aq. The preferred values are Oq ~ 0.3 and Aq ^ 0.7. 

1.3.6. Anthropic Considerations 

In models such as chaotic inflation, in which the observed big bang is just one of an infinite number 
of expanding regions in each of which the fundamental property takes a diff'erent value (Linde 1986, 1987, 
1988), and models in which a state vector is derived for the universe which is a superposition of terms with 
different values of the fundamental property (e.g. Hawking 1983, 1984; Coleman 1988), the probability of 
observing any particular values of the cosmological parameters is conditioned by the existence of observers in 
those "subuniverses" in which the parameters take these values (Efstathiou 1995; Vilenkin 1995; Weinberg 
1996; Martel, Shapiro, & Weinberg 1998). This probability is proportional to the fraction of matter which is 
destined to condense out of the background into mass concentrations large enough to form observers. Using 
this approach, Martel et al. (1998) calculated the relative likelihood of observing any given value of the 
cosmological constant A within the context of the flat CDM model normalized to COBE, and found that 
small but finite value of the cosmological constant, in the range suggested by observations, are favored over 
the value A = 0. Garriga, Tanaka, & Vilenkin (1998) have performed a similar analysis, but applied the the 
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density parameter, and found that intermediate values of fio are more likely to be observed than values near 
or near 1 . Anthropic arguments do not favor the values that the parameters take in the Standard Model, 
= 1 and Xq = 0. 

1.3.7. Alternatives to the Standard Model 

The problems listed above strongly argue against the Standard Model, and forces us to consider alter- 
natives. The age problem, large-scale structure problem, and the baryon catastrophe can all be solved by 
considering Open CDM, or OCDM models, in which the density parameter flQ < 1. However, such models 
do not satisfy the flatness requirement of inflation. CDM models with a nonzero cosmological constant Aq 
equal to 1 — Hq, known as ACDM models, satisfy this flatness requirement, and the addition of the cosmo- 
logical constant improves the age and large-scale structure problems, while providing a better agreement to 
the supernovae data. Recently, several authors have shown that it is possible to reconcile the inflationary 
scenario with an open universe, thus eliminating the flatness requirement (Ratra & Peebles 1994; Bucher, 
Goldhaber, & Turok 1995; Yamamoto, Sasaki, & Tanaka 1995; Linde 1995; Linde & Mezhlumian 1995). This 
not only supports open, matter- dominated models, but also allows for the possibility of an open universe 
with Ao 7^ and + Ao < 1. 

The large-scale structure problem can also be solved by introducing a "tilt" in the primordial power 
spectrum. In this Tilted CDM, or TCDM model, the primordial power spectrum P{k) at large scales does 
no have the Harrison-Zel'dovich form P(fc) oc fc, but instead varies as P{k) oc fc", where the primordial 
exponent n can differ from unity. The universe might also contain a mixture of two different forms of dark 
matter, one cold and one hot, a model known as CHDM. Finally, the universe might contain a smooth 
component whose pressure p and density p are related by an equation of state p = wp, a concept known 
as "quintessence" (Caldwell, Dave, & Steinhardt 1998: sec also Fry 1985; Charlton & Turner 1987; Silveira 
& Waga 1994; Martel 1995; Martel & Shapiro 1998). The cosmological constant is a particular form of 
quintessence, corresponding to w = —1; other forms have been suggested, such as domain walls, textures, or 
strings. 

The Standard Model had no free parameters. The values of Clo and Aq were fixed at 1 and 0, respectively. 

The value of h had to be close to 0.5 in order to avoid conflicts with ages of globular clusters, and the 
primordial power spectrum was assumed to be a CDM spectrum with no tilt. With the emergence of 
alternative models, there are now many free parameters. The density parameters fio no longer has to be 
unity. The cosmological constant Aq can be nonzero and, if open inflation is correct, a nonzero Aq does not 
have to be equal to 1 — Oq. The Hubble constant can vary over a certain range without conflicting with 
observations, and the slope n of the primordial power spectrum does not have to be 1. In models such as 
CHDM and quintessence models, there are additional parameters: the contribution of each component to 
the density parameter, and in the case of quintessence models, the coefHcient w appearing in the equation 
of state. 

During the 1980's, N-body simulations have played a central role in establishing the Standard Model, 
and then went more or less into hibernation, as efforts were invested into adding more physics to the original 
algorithms (hydrodynamics in particular). The emergence of alternative cosmological models has lead to a 
renewal of interest in N-body simulations. Such numerical simulations are essential for testing cosmological 
models against observations. Furthermore, they are useful from a theoretical viewpoint, since they can reveal 
how each cosmological parameter affects the process of large-scale structure formation. 
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1.4. The Need for a Database 

Numerical methods such as the Particle-Mesh algorithm (PM) and the Particle-Particle/Particle-Mesh 
algorithm (P'^M) are well documented. Details of the algorithms can be found in textbooks (e.g. Hockney 
& Eastwood 1981) and papers (e.g. Efstathiou et al. 1985). Hence, any researcher can easily access all the 
information and knowledge necessary to develop such algorithms. However, the effort required to develop, 
test, and optimize a PM or P^M algorithm from scratch can be quite substantial, and can be regarded 
as a waste of effort, since it essentially amounts to "reinventing the wheel." Also, performing simulations 
with large number of particles can demand a substantial investment in resources such as computer time, 
which is also wasteful if these simulations, or similar ones, have already been performed by other researchers. 
Consequently, it is a common practice among researchers to share either their programs or the results of 
their simulations. 

Klypin & Holtzman (1997) have combined into a single package their version of the PM algorithm and 
programs for generating initial conditions and analyzing the results. This package has been made available 
to the astronomical community, and can can be downloaded from a world-wide-web site. This allows other 
researchers interested in performing cosmological numerical simulations to "get started" immediately, without 
having to develop and test any computer program. However, installing and running these programs might 
pose some difficulties depending upon the kind of computer resources available to the user. 

We use a different approach. Instead of making our programs available to the astronomical community 
(something we might do eventually), it is the results of the simulations themselves that we are making 
available. We performed a very large number of numerical simulations, a total of 160, for 68 different 
cosmological models. This constitutes by far the largest database of cosmological simulations ever assembled, 
and it is still growing as more simulations are being performed. We are making this database available to 
the astronomical community (see §3.4 below). This approach is complementary to the one used by Klypin 
& Holtzman. By providing the results of the simulations, we eliminate the need for researchers to perform 
themselves these simulations, and the same simulations can be used by many different researchers. However, 
someone might be interested in simulating a cosmological model which is not included in the database, in 
which case the algorithm of Klypin & Holtzman can be used. Alternatively, we can, upon request, perform 
additional simulations and include them in the database. An interesting question is whether the results of 
simulations from the database can be analyzed using Klypin & Holtzman programs. In principle, this should 
be possible. The output files in the database are not written in the same format as the ones produced by 
Klypin & Holtzman's PM code, but it is fairly trivial to write a program that "translate" files from one 
format to another. 

There is an important difference that must be pointed out. The program of Klypin & Holtzman is 
based on the PM algorithm, while the simulations in the database were performed using a P'^M algorithm. 
For a same number of particles, the P^M algorithm has a length resolution superior to the one of the PM 
algorithm by a factor of order 6 (depending upon the particular choice of smoothing length). However, since 
the PM algorithm is significantly faster than the P'^M algorithm, it is possible to make up for the lack of 
resolution of the PM code by simply using more particles. We used 64'^ particles in all simulations]^ These 
have the same length resolutions as PM calculations with ~ 384'^ particles (such as the ones performed by 
Gross et al. [1998]). 

As mentioned above, there are numerous alternatives to the Standard Model. In this paper we consider 



^We intend to add simulations with 128'^ particles to the database in a near future. 
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CDM models in which the only components are ordinary matter (dark and baryonic) and possibly a nonzero 
cosmological constant (thus excluding CHDM models, and generic quintessence models). We consider the 
three cases f^o = 1, Aq = (the Einstein-de Sitter model), < 1, Aq = 0, and fio + Aq = 1- We also allow 
the primordial power spectrum to have a tilt. These models are usually referred to as Tilted CDM (TCDM), 
Tilted, Open CDM (TOCDM), and Tilted, Lambda CDM (TACDM). 



2. THE NUMERICAL SIMULATIONS 
2.1. The Algorithm 

All simulations presented in this paper were done using the P^M algorithm (Hockney & Eastwood 1981; 
Efstathiou ct al. 1985). The computational volume is a cubic box of comoving size I/box and comoving volume 
14ox = ibox with triply periodic boundary conditions, expanding with Hubble flow. The matter distribution 
inside the computational volume is represented by N equal-mass particles. The forces on particles are 
computed by solving Poisson's equation on a cubic grid using a Fast Fourier Transform method. The 
resulting force field represents the Newtonian interaction between particles down to a separation of a few 
mesh spacings. At shorter distances the computed force is significantly smaller than the physical force. 
To increase the dynamical range of the code, the force at short distance is corrected by direct summation 
over pairs of particles separated by less than some cutoff distance re- With the addition of this so-called 
short-range correction, the code accurately reproduces the Newtonian interaction down to the softening 
length e, which is a fraction of the grid spacing. The system is evolved forward in time using a second order 
Runge-Kutta time-integration scheme with a variable time step. 

Our particular version of the P^M algorithm uses supercomoving variables (Martel & Shapiro 1998; see 
also Shandarin 1980). In these variables, the position r, peculiar velocity v, time i, density p, and peculiar 
gravitational potential (j) are related to their Eulerian counterparts by 



avi* 



(5) 
(6) 



p* 
a^tl 



(9) 



where 



P* = Po= , (10) 

2 



In these equations, a{t) is the Priedmann- Robertson- Walker scale factor, ao is its present value, and r*, is 
a free parameter whose value is chosen according to the characteristic length scale of the problem. These 
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variables are similar to the more standard comoving variables in many respects. In particular, equations 
(5) and (8) imply that a volume expanding with Hubble flow remains fixed in supercomoving variables, and 
that the mean density inside that volume remains constant. The main difference is in the change of time 
variable given by equation (7). In supercomoving coordinates, the time i is negative, and equal to — oo at 
the big bang. In an Einstcin-dc Sitter model, f = — 1 at present. This change of time variable has the virtue 
of eliminating the cosmological drag term in the momentum equation. 

In all simulations, we set = Lbox/cio- Equation (5) then implies that the box size in supercomoving 

variables is unity at all times. 

The time-evolution of the scale factor a{t) is governed by the Friedmann equation. For universes 
composed of ordinary, nonrelativistic matter and a nonzero cosmological constant Aq, the Friedmann equation 
takes the form 



1 da^ ^ 



(l-f2o-Ao){ +"°(^^ 



(12) 



In supercomoving variables, there is a precise normalization for the scale factor, which depends upon the 
particular cosmological model. For the models considered in this paper, the solution of the Friedmann 
equation and the present value of the scale factor are the following: 

(a) Einstein-de Sitter model {Qq = 1, Aq = 0) 

a = i~'^ , ao = 1 . (13) 

(b) Open models (Qq < 1, Aq = 0) 

a = {P - l)-i , ao = (1 - no)/no . (14) 

(c) Flat models with nonzero cosmological constant (Oq + Aq = 1) 

1 d,y ( ^ \ ^^"^ 

Notice that the solutions for a{t ) do not depend explicitly upon the cosmological parameters, which are 
absorbed in the definition of ao. Hence, for all models included in the database, there are only 3 different 
solutions of the Friedmann equation. This is one of the most useful properties of supercomoving variables. 
For simplicity, we shall drop the tilde notation for supercomoving variables in the remainder of this paper, 
except in §3.3. 



2.2. The Power Spectrum 

For all simulations presented in this paper, we use the Cold Dark Matter (CDM) power spectrum of 
Bardeen et al. (1986), with the normalization of Bunn & White (1997). The power spectrum at redshift z 
is given by 

P(fc, z) = 27r' i^—j S%C-\z, 0)fc"T^DM(fc) , (16) 
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where c is the speed of light, C{z,0) = S+{0)/S+{z) is the Hncar growth factor between redshift z and the 
present, and S+ is the Unear growing mode (see eqs. [35]-[39] below), n is the tilt, and Tcdm is the transfer 
function, given by 

TcdmW = ^"^^+3^^^^'^ [1 + 3.89g + {IG.lqf + {5A6qf + (6.71g)^] "^/^ (17) 
(Bardeen et al. 1986), with q is defined by 

5 = f77^V"'^'("o'^')"'©2.7, (18) 



VMpc 

a = a-"Bo/no^-(nBo/no)^ ^ ^^g^ 

ai = (46.9r!o/i')°-^^"[l + (32.ir!o/i2)-°-^32] , (20) 
as = (12.0J7o/i')°-^'^[l + (45.017o/i')-°-H (21) 

(Hu & Sugiyama 1996, eqs. [D-28] and [E-12]), where 62.7 is the temperature of the cosmic microwave 
background in units of 2.7K. and Sh is the density perturbation at horizon crossing (Liddle & Lyth 1993). 
Fits for 6h are given by Bunn & White (1997), as follows, 

in5, _ ri.95fi,7""'«'"""-""^e-(^+"i4*i^), Ao = 0; 

" I 1.94f]-0-785-0 051„nog-(0.95n+0.169n^) ^ Aq = 1 - f^o; ^ ' 

where n = n — l. 



2.3. Setting up Initial Conditions 

We assume that the initial fluctuations originate from a Gaussian random process. The initial density 
contrast can then be represented as a superposition of plane waves with random phases, and amplitudes 
related to the power spectrum P{k), where k is the wavenumber, and k = |k|. In an infinite universe, all values 
of k arc allowed. The power spectrum is therefore continuous, and the number of modes (that is, plane waves) 
present in the initial density contrast is infinite. The simulations, however, are performed inside a finite 
comoving cubic volume Vbox = -^box with periodic boundary conditions. This periodicity implies that only 
modes with wavenumbers k = (fc^, ky, kz) = {nx,ny,nz)ko, where n^, riy, Uz are integers and ko = 27r/Lbox 
is the fundamental wavenumber, can be present in the simulated initial conditions. Furthermore, since 
the initial conditions are represented by particles, the components kx, ky, kz of the wavenumber cannot 
exceed the nyquist frequency ftnyq = iV^/^fco/2, where N is the number of particles in the computational 
volume, and 7V^/^ is the number of particles along one dimension. Modes with higher wavenumber cannot 
be represented because of imdcrsampling. Hence we are faced with the task of representing continuous initial 
conditions using a discrete sample of plane waves. This key aspect of any numerical cosmological simulation 
is, surprisingly, seldom discussed in the literature. Here we present a detailed description. 

In a periodic universe with comoving cubic volume Vbox = -^box' density contrast 6 can be decom- 
posed into a sum of plane waves. 

Sir) =Y,St''e-^''-'' , (23) 

k 
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where r is the comoving, or supercomoving, position, and J^"*"^ is the ampUtude of the mode with wavcmimbcr 
k. The superscript "disc" stands for "discrete." The real universe is of course not periodic, in which case all 
values of k are allowed. To convert equation (23) from the discrete limit to the continuous limit, consider 
first any function /(k) that is summed over all possible values of k. In the discrete limit, we have 

k allV.E. '^0 allV.E. ^^^> all V.E. "^V.E. 

where "V.E." represents a volume element va k-spacc, which is a cube of volume fcg = Vbox/(27r)'^ centered at 
k. Assuming that the function / does not vary significantly over one volume element, we can pull it inside 
the integral, 

E/k'^^«7CT E / f^'^'^"^- (25) 



(2-) allV.E. 



V.E. 



Of course, integrating over the volume element, and then summing over all volume elements, is effectively 
like integrating over all k-space, so equation (25) reduces to 



(26) 



where the superscript "cont" stands for "continuous." The continuous and discrete functions are related by 

/r* = 7^/k"^ (27) 



Using these formulae, we can rewrite equation (23) as 

S{r) = J rf3fc^£°"*e-*''•^ (28) 

where 



rconl ^ox cdisc /'on^ 



To find the relationships between J^''*^, and the power spectrum, consider the rms density fluctu- 

ation ax at some particular scale x. This quantity is given by 



= (^ y d'k\StrWikx) . (30) 

where is a window function. We present the derivation of this result in Appendix A. In the continuous 
limit, is related to the power spectrum, by 

<jI = ^ j £kP{k)W{kx) (31) 

(see, e.g. Bunn & White [1997], eqs. [22] and [23]). By combining equations (29), (30), and (31), we get 

P{k) = V^o.lSt'l^ = ^l-^rf • (32) 
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Both P{k) and have dimensions of a volume while S^^^"^ is dimensionless. Notice that the form of these 
expressions depends upon the actual definition of the Fourier Transform, which tends to vary among authors. 

To set up initial conditions, we lay down the particles on a cubic lattice, and displace each particle by 
an amount Ar given by 



where r is the unperturbed position, 6^^^'^ = \S^^^'^\e~'^'^^ is a complex number with amplitude \S^^^^\ = 
[P(fc)/Vbox]^^^ and phase 0k chosen randomly between and 27r with uniform probability, and the sum 
extends over all modes included in the initial conditions (sec §2.4). As Efstathiou ct al. (1985) point 
out, assuming random phases would be sufficient to ensure that the initial conditions are Gaussian, in the 
continuous limit (that is, in an infinite universe). However, this assumption is insufficient in the discrete 
limit (that is, in a finite universe with periodic boundary conditions). To ensure the Gaussianity of the initial 
conditions, it is necessary, and sufficient, to include the Gaussian factor Gk, a random number chosen from 
a Gaussian distribution with mean and dispersion 1. This guarantees the initial conditions are Gaussian, 
even though there might be a lack of resolution at some scales; Gk does not change the spectral amplitude 
of the fluctuations. 

To compute the initial peculiar velocity fleld, we assume that the initial time of the calculation is early 

enough for the perturbation to be in the linear regime, but late enough so that the linear decaying mode 
can be neglected. The initial peculiar velocity of the particles are then related to their displacements by 



where Zi is the initial redshift of the simulations, Ar is computed using equation (33), and i5+ is the linear 
growing mode of the perturbation, which depends upon the cosmological model. For the Einstein-de Sitter 
model {fto = 1, Aq =0), the growing mode is 




(33) 



k 




(34) 



6+{z) = {l + z) 



-1 



(35) 



For open models (l^o < 1, Aq = 0), the growing mode is 




(36) 



(Peebles 1980), where 



X = {^0^ - + z) 



1 



(37) 



Finally, for flat models with a cosmological constant {flo + Ao 



1), the growing mode is given by 




(38) 



(Martel 1991b), where 



(39) 
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2.4. 



The Simulations 



We set the comoving length of the computational volume ibox equal to 128 Mpc (present length units). 
The total mass of the system is M^ys = SH^QoLI^^/SttG = 5.821 x lO^'^f^o/i^MQ. We use = 64^ = 262, 144 
particles of mass Mpart = Msys/N = 2.220 X lO^^f^o/i^M©. We solve Poisson's equation on a 128"^ grid. In all 
simulations, e and re were set equal to 0.3 and 2.7 mesh spacings, respectively. This corresponds, in physical 
units, to a comoving softening length e = 300 kpc. This is a reasonable value for gravity-only cosmological 
simulations. At smaller scalcis. hydrodynamical effects become important and cannot be ignored. The 
dynamical range in length of the algorithm is Lbox/e = 467. 

The ratio of the nyquist wavenumber fcnyq to the fundamental wavenumber fcg is N^^^ /2 — 32. Hence 
each component /cj, i = x, y, z, of the wavenumber can take 65 values; ki = riiko, with —32 <ni< 32. The 
initial conditions can therefore represent 65^ = 274 625 modes. However, the reality condition requires that 
the amplitudes of modes with equal and opposite wavenumbers are related by 5^^^^ = {S"^^)* [in order for 
(5(r) to be real]. Furthermore, we exclude modes with |k| = {k^^ + ky + k'^Y^^ > fcnyq. This reduces the actual 
number of modes represented in the initial conditions to 68 532. 

All simulations start at an initial redshift Zj = 24. The algorithm produces "dumps" (snapshots of the 
system) at numerous intermediate redshifts, up to the present. These redshifts where chosen by imposing 
that the dumps are equally space in conformal time ij, defined by dij = aodt/a{t). We set the difference Ary 
between consecutive dumps equal to ibox/ c. Thus, if t and t' are the times corresponding to 2 consecutive 
dumps, they are related by 



This particular choice results in most dumps being concentrated near the present. Typically, about half of 

the dumps are between redshifts z = 1 and z = 0. Since the relationship between time and redshift, z{t), 
is model-dependent, the redshifts where dumps are made depend upon the cosmological parameters flo, Aq, 
and Ho (but not as)- Every simulation also produces a dump a,t z = Zi = 24, and one at 2: = 0. The number 
of dumps per simulation varies between 44 and 128. 



The power spectrum described in §2.2 is characterized by 6 independent parameters: (1) the density 
parameter f^o, (2) the contribution fieo of the baryonic matter to the density parameter, (3) the cosmological 
constant Aq, (4) the Hubble constant Ho, (5) the temperature Tcmb of the Cosmic Microwave Background, 
and (6) the tilt n of the power spectrum. In order to keep the size of the parameter space at a manageable 
level, we set Tcmb = 2.7 K and fJeo = 0.015/i~^, thus reducing the dimensionality of the parameter space to 
4. Also, the normalization of the power spectrum is often described in terms of the rms density fluctuation 
as at a scale of 8/i~^Mpc. The value of as is a function of the 6 aforementioned parameters. We invert 
this relation, treating as as an independent parameter, and the tilt n as a dependent one. The independent 
parameters in the database are therefore Qq, Aq, Hq, and as- For each model, we performed up to 3 different 
simulations, with different realizations of the initial conditions (this amounts to choosing a different set of 
random numbers for the phases (^k of the complex numbers 6^^^^'^, and the Gaussian factors Gu). 




(40) 



3. 



THE DATABASE 



3.1. The Cosmological Models 
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An important question was to decide which models should be included in the database. Our goal here 
is not to find the "ultimate model," which provides the best match to current observations. This would 
defeat the purpose of having a database, and furthermore, as former supporters of the Standard Model can 
appreciate, the "best model" can eventually be proven incorrect by new observations. Our intention is to 
provide an adequate coverage of the parameter phase-space. However, we do not want to invest much effort 
into simulating models that are considered "unlikely," because some of the parameters have extreme values. 
With this in mind, we performed 160 simulations, which provide a broad coverage of the parameter phase- 
space, but we favored "likely" regions of the parameter phase-space over "unlikely" ones, by performing more 
simulations in these regions. For instance, we consider models with Hubble constant varying in the range 
i?o = 50 — 85 kms~^Mpc~^, but 139 of the calculations (87%) have a Hubble constant in the more plausible 
range Hq = 65 — 75kms~^Mpc~^. 

The value of the parameters are given in Table 1 for the entire database (with Hq in units of km s~ ^ Mpc~ ^ ) . 
The first 4 columns contain the values of the 4 independent parameters Qq, Aq, Hq, and erg. The dependent 
parameter n is in the fifth column. The sixth and seventh columns contain the number of dumps per sim- 
ulation and the codes of the simulations respectively (see §3.2). The parameter phase-space coverage of the 
database is illustrated in Figure 1. The top left panel shows a projection of the 4-dimensional parameter 
phase-space onto the flo — Aq plane. The dots indicate the cases for which there are simulations in the 
database. The number next to each dot indicates the number of simulations for that particular combination 
of and Aq. This panel includes all simulations in the database. The top right panel shows the same 
projection, but for a subset of the simulations, all simulations with Hq = 65kms~^Mpc~^. The remaining 
4 panels show different projections and different subsets. As we see, the coverage of the parameter phase 
space is quite dense. The biggest "hole" is seen in the flo — Aq projection (top panels). There are currently 
no simulations for open models with a nonzero cosmological constant (Aq 7^ and fio + Aq < 1) in the 
database. As we pointed out in §1.3, these models are certainly worth considering, and we intend to include 
such models in the database in the near future. 

Several interesting quantities can be computed directly from the parameters. One of them is the age of 
the universe. For Xq ^ models, to is given by 



Aoa;'^ -h (1 — — Xq)x + flox 




(see, e.g., Martel 1990). It is of course independent of erg. Table 2 gives the ages is Gigayears for the various 
models included in the database. 

Another interesting quantity is crf^^, the value of erg inferreded from observations of clusters of galaxies. 
Using the X-ray temperature distribution function of clusters, Viana & Liddle (1996) have produced an 
empirical formula for erg'"'*, 

£7^'"^ = Om~'^^^°^ , (42) 



where 



0.36 -I- 0.31f2o - 0.280^ , Aq = ; 
0.59 - 0.16^0 + 0.06^^ , Ao = 1 - Oo . 



<-'l^^0) - i ^ ^„ , r,nco2 x _1 r, (43) 



Table 3 gives the values of af^^ for the various models included in the database. Notice that these values 
do not always match the actual values of erg used for the calculations (fourth column of Table 1). 



3.2. Nomenclature 



Each cosmological model, that is, each combination of the four parameters CIq, Aq, Hq, and ug, is 
identified by a two-character code composed of an uppercase letter and a lowercase letter. For instance, the 
Einstein-de Sitter model with Hq = 65kms~-'^Mpc~^ and cts = 1-0 is identified by the code Xa. The letters 
were chosen for practical reasons, and the reader should not try to find some logic in these choices. Each 
simulation is identified by a 3-character code, composed of the two-character code for the model, plus a digit 
to identified the simulation. For instance, the three simulations for the Xa model are identified by the codes 
Xal, Xa2, and Xa3. The codes for the entire database are given in the last column of Table 1. Each simulation 
produces many output files, or dumps, which are snapshots of the system at various redshift. Each file is 
identified by a 7-character code, which consists of the 3-character code of the simulation, an underscore, and 
a 3-digit number which identifies the file. For instance, the first output file created by the simulation Xal is 
called Xal_001, and contains a snapshot of the system at the initial redshift Zi = 24. The next file created 
by that simulation is called Xal_002, and contains a snapshot a.t z = 22.079, and so on. The last file is called 
Xal_058, and contains a snapshot at the present {z = 0). The lists of redshifts where dumps are available 
can be obtained from the authors. There is one such list for every combination of flo, Aq, and Hq included 
in the database. 



3.3. Conversion to Physical Units 



The positions and velocities stored in the dumps are expressed in supercomoving variables. They can 
be converted to physical units using equations (5)-(ll). The positions r and velocities u = H{z)r -|- v in 
physical units are given by 



-t'boxf 
U = Lbox 



H{z)r _^i^^Ho{l + z)v 



1 + z 



(44) 
(45) 



In these expressions, we have reintroduced the tilde notation for the supercomoving variables. The expression 
of r is the same for all models, but the one for v is model dependent. After eliminating Ji{z) using 
equation (12) and ao using equations (13)-(15), we obtain the following expressions: 



(l+^)l/2f + 



(a) Einstein-de Sitter model (Oq = 1, Aq = 0) 

(b) Open models (fio < 1, Aq = 0) 

u = ifo-£'box (1 -h Oo2;)^/2f -h 



(l+z)v' 



no(l + -z)v 



2(1 - Oo)V2 

(c) Flat models with nonzero cosmological constant (f2o + Ao = 1) 

u = i^oLboxj [^^0(1 + + Ao] ^/^^ ^"^^^Y 



(46) 



(47) 



(48) 



3.4. Technical Considerations 



The database contains 160 simulations for 68 cosmological models. For each simulation, there is a 
dump at z = Zj = 24 and one at z = 0, plus numerous dumps at intermediate redshifts. There is a total 
of 11973 dumps in the database. Each dump contains 6A^ = 1572 864 numbers, the coordinates of the 
position and velocity for each particle in the simulation. These numbers are stored in single precision (32 
bits), although the simulations themselves were performed in double precision. Each dump is a binary 
file (IEEE 754 standard) of size 6.3 Mb, which contains first the x-coordinatcs of all particles, followed by 
the y-coordinates, the z-coordinates, the t;cc-coordinates, the ?;y-coordinates, and finally the ?;z-coordinates. 
Figure 2 shows a sample FORTRAN program that reads a file from the database. The size of the entire 
database is 75.3 Gb. The database currently resides on archival tapes at the High Performance Computing 
Facility, University of Texas, where the simulations were performed. 

Because of the size of the database, it would be impractical (if not impossible) to install it on a web 
site or an anonymous ftp site where it could be easily retrieved by the user. This might change in the 
future, but currently the only way to access the database is to contact the authors, preferably by E-mail, 
and send a list of the dumps requested. Then, the authors and the user can choose the best strategy for 
transferring file, according to the computer resources and needs of the user. Requests should be sent to 
databaseQgalileo . as . utexas . edu. 



The goal of this paper is to present the database, and describe its content. We supplement this descrip- 
tion by analyzing the final state of each simulation, which corresponds to the present. We focus on four 
particular aspects of the present large-scale structure: the rms density fluctuation, the two-point correlation 
function, the moments of the peculiar velocity field, and the properties of clusters. In this abbreviated ver- 
sion of the paper, we only present the analysis of the rms density fluctuation and the two-point correlation 
function. The full version of the paper can be obtained by contacting the authors. 



The present rms density fluctuation as at scale 8/i~^Mpc is treated as an independent parameter in 
our simulations. However, we do not set up the state of the system at present. Instead, we set up initial 
conditions at high redshift [zi = 24), and evolve the system numerically all the way to the present. We 
adjust the initial conditions in such a way that the density fluctuation at presents ends up being equal to 
the desired value of a^. To achieve this, we assume that the power spectrum evolves with time according 
to linear perturbation theory (hence the presence of the factor in equation [16]). Actually, we do not 
expect the actual value of erg to be precisely equal to the desired one, for several reasons. Let us designate 
by cTg""* the desired value of (Tg for each simulation, that is, the quantity appearing in the fourth column of 
Table 1. The superscript "cont" indicates that this the value in the real universe, where the wavenumber 
k varies continuously. We are representing the initial conditions using a finite number of discrete modes, 
which is clearly an approximation. We designate by crg'^'^ the value of ag resulting from this approximation. 
This value is given by 



4. 



ANALYSIS OF THE SIMULATIONS 



4.1. RMS Density Fluctuation 




(49) 
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where £ = 8h ^Mpc. Hence, cTg^^'' is computed by summing over all modes present in the initial conditions, 
ignoring the Gaussian factor, and then extrapolating to the present using linear perturbation theory. 

The introduction of the Gaussian factor in equation (33) further modifies the value of erg. We designate 
by o'l''"'''' the value of erg resulting from the presence of this factor, 

(ar-)2 = C\z,, 0) ^ Gl\6if\''Wike) . (50) 

k 



Finally, we designate by ct""™ the "numerical" value of as , which is the actual rms density fluctuation 
inside the computational volume at present, obtained from the numerical simulation. This value should 
differ from af"^^ from several reasons. First, the numerical algorithm has a finite accuracy, owing to the 
fact that the time step is finite and that the gravitational force is softened at short distances. Second, the 
evolution of a single mode would never follow precisely the exact solution when the system is represented by 
a finite number of particles. Third, and more importantly, equations (49) and (50) use linear perturbation 
theory to extrapolate from the initial conditions to the present, but this can only be approximate, as mode 
coupling introduces nonlinear effects at small scale. 

We investigated the importance of these various effects by computing the various values of as for 
all simulations. The value of cr™"* is imposed. The vahics of a^^^^ and erf''"'''' can be computed directly 
using equations (49) and (50) (these values are provided automatically by the code that generates initial 
conditions). To evaluate erg"™, we used a direct, somewhat brute-force approach. For each simulation, we 
located one million (!) spheres of radius 8/i~^Mpc at random locations inside the computational volume at 
present, and computed the density contrast Jgph inside each sphere, using 

4ph = ^^HP^, (51) 

-/Vsph 

where A''sph is the number of particles inside the sphere, and Nsph is the "mean number," given by 

A^sph = (52) 

(notice that iVgph is not an integer). The value of ag"™ is then given by 



(a8-™)2 = 10-6 ^ («^ (53) 

all spheres 



We plot these various values of as against each others in Figure 3. The top-left panel shows the effect of 
discreteness. All dots are located below the dashed line, indicating that crf'^'^ < cr™"*. This is caused not as 
much by the discreteness itself as by the fact that modes outside the range k^ ^ k ^ ^nyq are missing in the 
initial conditions. Still, the effect is very small, 7% in the worst case (which happens to be model Uc). As we 
see on the top right panel, the effect of introducing of the Gaussian factor Gk is quite important, and causes 
a spread of order 10% in the value of erg. This is primarily an effect of undersampling. The modes with 
wavenumbers comparable to fco are very few, but contribute significantly to erg. This is consequence of the 
fact that our computational volume is actually too small to constitute a "fair sample" of the universe. With 
a larger volume, there would still be very few modes with wavenumbers of order ko, but these modes would 
be farther from the peak of the power spectrum, and therefore would contribute less to erg. The bottom left 
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panel shows the effect of actually performing the simulation. The values of frf'^"'''' and ct™™ are comparable 
for small cts, but differ at large cts where nonlinear effects become important. This panel shows that (1) 
the onset of nonlinearity occurs at erg ~ 0.6, (2) The effect of nonlinearity on the value of CTg is small, of 
order 10%, and (3) the effect can go either way, it can either increase or decrease the value of erg, and the 
occurrences of these two cases are comparable. 

The bottom right panel shows the combined effect of discreteness, Gaussian factor, and nonlinearity. 
The spread is quite large. The value Cg"™ that comes out of the simulation can differ by as much as 20% 
from the value cr§°'^* we intended to obtain, that is, the values given in the fourth column of Table 1. The 
reader should be aware of this fact when selecting a particular model from the database. Consequently, we 
listed in Table 4 the values of a-?"™ for all simulations. 



4.2. 2-point Correlation Function 

The distribution of galaxies in the universe can be described statistically using N-point correlation 
functions. The first and most important of these functions is the 2-point correlation function ^(r), which 
measures the excess probability of finding two galaxies separated by a distance r. The 2-point correlation 
function can be estimated from the particle distributions by locating spherical shells around particles and 
counting the number of particles inside these shells. If N{ri, is the number of particles inside a spherical 
shell of inner radius ri and outer radius r2 centered on a given particles, and {N{ri,r2)) is the average of 
N{ri,r2) over all particles, then by definition 

{N{ri,r2)) = ^^(""z ~ '^?)'" + 47r„ H ^{rydr . (54) 



where n is the number density of particles. The first term in equation (54) gives the correct answer in the 
case of a uniform distribution. The second term represents the effect of the correlation. We can solve this 
equation for ^ (see Martel 1991a). After some algebra, we get 



r^{x) In 10 da; 



{N{r{x),n)) r^{x) 



Aim 



(55) 



where x = logj^p To compute ^, we first evaluate {N{r{x),ri)) by computing the spacing between all pairs 
of particles, and counting pairs in bins equally spaced in x. Wc then compute the derivative in equation (55) 
using a standard five-point finite difference operator. We computed ^ (r) at present for all simulations in the 
database. The results are plotted in Figures 4-7. For models with more than one simulation (most have 
thrc;e), wo averaged the curves. They were actually so similar, in all cases, that error bars in Figures 4—7 
would be too small to be seen. This shows that the 2-point correlation function depends essentially upon 
the cosmological model, with very little dependence upon the particular realization of the initial conditions. 
Each panel in Figures 4-7 corresponds to a particular combination of Oqj Aq, and Hq, with different values of 
cTg represented by different curves. The dashed lines show the observed galaxy 2-point correlation function, 

-1.77 



(Peebles 1993, eq. [7.32]). Models with erg = 0.3 fail to reproduce the observed correlation function on three 
counts. The slope of ^(r) is too shallow, the amplitude is too small, and there tend to be a "kink" in the 
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correlation function at separations of order 1 — 3 Mpc. As a% increases, the kink goes away, and the amphtude 
and slope increase. Models with a% — 0.8 provide the best fit to the observations, almost independently of 
the values of the other parameters! For larger values of erg, the slope and amplitude arc too large. 

All curves have a shoulder at small separations, where the slope drops significantly. Martel (1991a) 
argues that this is a consequence of the softening of the force at small distances. The effect of this softening 
is to "take" pairs of particles that would have a separation r less than the softening length c = 300 kpc in 
the absence of softening, and transfer them to separations r ^ e, resulting in a flattening of the correlation 
function. To illustrate this, we indicated in all panels of Figures 4-7 the location of the softening length e 
by a thick line. This line is located right in the middle of the "shoulder" for all curves. With higher force 
resolution, pairs that are now located at separations r ^ e would be located instead at separations r <^ e, 
and the fit to the observed slope would be improved. 

5. SUMMARY AND PROSPECTS 

Using a P^M algorithm with 64^ particles, we have performed 160 cosmological simulations, for 68 
cosmological models. This constitutes the largest database of cosmological simulations ever assembled. We 
covered a four-dimensional parameter phase space by varying the density parameter fio, the cosmological 
constant Aq, the Hubble constant Hq, and the rms density fluctuation erg. We are making this database 
available to the astronomical community. We also performed a limited analysis of the simulations. Our 
results are the following: 

(1) The present rms density fluctuation erg"™ differs from the one expected by linearly extrapolating the 
initial power spectrum to the present, because of the combined effects of having a finite number of modes in 
the initial conditions, introducing a Gaussian factor in the initial conditions, and having nonlinear coupling 
between modes. The first effect is negligible. The second effect is of order 10% or less, but this probably 
depends upon the size of the computational volume [the one used for the simulations, (128 Mpc)'', is a little 
too small to constitute a fair sample of the universe]. The effect of nonlinearity is negligible for ag < 0.6, 
and of order 10% or less for larger erg. This can go either way: nonlinearities can either increase of decrease 
(Tg relative to what linear theory predicts. 

(2) The observed two-point correlation function ^(r) is well reproduced by models with erg = 0.8, nearly 

independently of the values of the parameters CIq, Xq, and Hq. For models with (Tg < 0.8, the correlation 
function is too small, its slope is too shallow, and it has a kink at separations r = 1 — 3 Mpc. For models 
with cTg > 0.8, the correlation function is too large and its slope is too steep. 

(3) At small separations, r < 1 Mpc, the velocity moments satisfy the relations |Vr| « Hor and Vpp ~ 

indicating that small clusters have reached virial equilibrium. At larger separations, |Va| increases 
above the Hubble velocity, indicating that clusters are accreting matter from the field. The velocity moments 
depend essentially upon Oq and erg, and not Aq and Hq. The pairwise particle velocity dispersions are much 
larger than the observed pairwise galaxy velocity dispersion, except for models with Oq = 0.2 and erg < 0.4. 
But if the velocity dispersion of galaxies is biased relative to the velocity dispersion of dark matter, then 
models with larger values of flo or cg can be reconciled with observations. 

(4) The multiplicity functions are decreasing for small values of for models with erg ~ 0.3. At larger values 
of erg, the multiplicity functions have a horizontal plateau, whose length increases with erg. For models with 
erg > 0.9, the multiplicity functions have a U shape which results from the merging of intermediate-size 
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clusters. For all models, clusters have densities in the range lOOpo — lOOOpo- A simple analytical model 
suggest that clusters have a density p ^ 178p when they reach virial equilibrium. Our results suggest that 
many clusters have reached that equilibrium in the past, when p was larger than po (this could be checked 
by performing a cluster analysis on earlier dumps). The spin parameters A are in the range 0.008 — 0.2, with 
the median near 0.05, and the distributions of elongations favors prolate shapes (e2 > ei) over oblate shapes 
(ei < 62). These results indicate the absence of rotationally supported disks in these simulations. 

The database is growing. We are currently adding new simulations to the original 160 simulations 
described in this paper. There are at least seven different motivations for performing additional simulations. 

(1) Additional Simulations for the same Models: For the sake of providing a good coverage of the parameter 
phase-space, we have limited the number of simulations per model to 3 or less. We can perform additional 
simulations for models already included in the database, if there is a need for doing so. This could be the 
case if, for some reason, a particular model (that is, a particular combination of fio, Aq, i?Oj ^^nd cg) becomes 
particularly interesting, and deserves more scrutiny. Also, having more simulations per model has the virtue 
of improving the statistics. For instance, the size of the error bars in Figures 8-12 would be reduced if we 
had more than 3 simulations per model. Finally, for gravitational lensing simulations, it is necessary to 
combine dumps generated by different simulations, and having more simulations can be desirable (see, e.g. 
Premadi, Martel, & Matzner 1998). 

(2) Different Box Size: All simulations included in the database were performed using a computational 
box of size ibox = 128 Mpc. As described in §2.4, the softening length is comparable to the scale where 
nongravitational effects become important. Hence, there is no reason to consider smaller box sizes, unless 
we want to use fewer particles. There are, however, reasons for considering larger boxes. As we pointed out 
in §4.1, a box of size 128 Mpc is too small to constitute a "fair sample" of the universe. Using boxes of size 
256 Mpc or even 512 Mpc would certainly provide a better, "fairer" description of the large-scale structure, 
even with the same number of particles. 

(3) Larger Number of Particles: The is no point increasing the number of particles as long as we keep the 
box size at 128 Mpc, since the resolution would be increased at scales where nongravitational effects are 
important. However, if larger boxes are used, the number of particles can be increased accordingly in order 
to maintain the resolution of the algorithm at small scale. If we continue to adopt 300 kpc as the resolution 
scale of the algorithm, simulations in (256 Mpc)^ and (512 Mpc)^ boxes could be performed with 128^ and 
256'^ particles, respectively. 

(4) New Background Models: The 4-dimensional parameter phase-space considered in this paper is quite 
large, and the set of 68 cosmological models included in the database covers a small fraction of it. There 
are several "holes" in the projections shown in Figure 1. In particular, there are no simulations for open 
models with a nonzero cosmological constant (Aq 7^ 0, fio + ^^o < !)• Simulations for additional background 
models could be added to the database, either to provide a better coverage of the parameter phase-space, 
or because there is a particular model we are interested in, "we" designating either the authors, or other 
researchers sending us a special request. Actually, the original database contained only 151 simulations 
for 65 cosmological models. Following a special request by Hamana (1998), we added 9 simulations to the 
database, for 3 new models: Ea, Pa, and Xg. 

(5) Additional Parameters: The current database covers a 4-parameter phase space, because we held the 
CMB temperature Tcmb and the baryon density parameter I^bo at values of 2.7 and 0.015/i~^, respectively. 
The CMB temperature is known so accurately that treating it as a variable parameter would be pointless. 
This is not the case for the baryon density parameter. According to primordial nucleosynthesis, the quantity 
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fieo^^ has an allowed range from 0.01 to 0.026 (Krauss & Kernan 1995; Copi et al. 1995; Krauss 1998). 
Furthermore, X-ray observations of clusters of galaxies suggest that the ratio of gas mass to dark matter mass 
in these clusters exceeds the mean value in the universe, a phenomenon known as "the baryon catastrophe" 
(Briel et al. 1992; White et al. 1993; Martel et al. 1994). There is at present no definitive explanation for 
this phenomenon, but one possible explanation is that primordial nucleosynthesis is somehow incorrect, and 
predicts a value of Qbo which is too small. 

(6) Different Components: All simulations in the database used a CDM power spectrum as initial conditions. 
There are, however, several other models that constitute interesting alternatives to the CDM model, which 
could be added to the database. One of them is the Hot Dark Matter model (HDM), though this model 
has fallen out of favor in recent years, duo to its inability to form galaxies inside deep voids such as Bootes. 
A more interesting alternative is the mixed Cold + Hot Dark Matter model (CHDM), which contains both 
a cold dark matter component and a massive neutrino component. This model introduces one additional 
parameter, the contribution fi^o of the neutrinos to the mean energy density of the universe. 

(7) Different Cosmologies: The cosmological models included in the database contain only non-relativistic 
matter and a nonzero cosmological constant. It would be very interesting to consider models with other 
components. Possible candidates include domain walls, cosmic strings, or relativistic particles (Fry 1985; 
Charlton & Turner 1987; Silveira & Waga 1994; Martel 1995; Martel & Shapiro 1998). Recently, these various 
candidates have been combined into a single concept called "quintessence" (Caldwell, Dave, & Steinhardt 
1998). The effects of these various components is twofold: First, the presence of these components modifies 
the expansion rate of the universe and the growth rate of density perturbations, thus changing the history 
of large-scale structure formation. Second, they might affect the shape and normalization of the primordial 
power spectrum, in ways that remain to be determined (none of these models were considered by Bunn & 
White [1997]). 
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edge the support of NASA Grants NAG5-2785, NAG5-7363, and NAG5-7821, NSF Grants PHY93 10083, 
PHY98 00725 and ASC 9504046, the University of Texas High Performance Computing Facility through the 
office of the vice president for research. HM acknowledges the support of a fellowship provided by the Texas 
Institute for Computational and Applied Mathematics. 



A. Calculation of C7^'"'= 
The mass inside a sphere centered at ro is given by 



M(ro) 



Pcom(l + S)(fr = Pc 



sph(ro) 



sph 



Jsph(ro) 1, 



isc — «k-r 

k e 



(Al) 



where pcom is the average comoving density, Vaph is the volume of the sphere, and the integral is computed 
over that volume. The relative mass excess in the sphere is given by 

AM, 



M ^ Vsph Jsph(ro) 



iph 7sph(ro) 1, 



(A2) 



We introduce the following change of variables. 



r = ro -I- y . 



(A3) 
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In y-space, the sphere is now located at the origin, and equation (A3) becomes 



jdisCg~ikrOg-ik-y 



(A4) 



We now square this expression, and get 



/AMV 



(ro) 



167r2x6 



Jsph(O) 



disc — ik-ro _— ik y 



Jsph(O) w 



(A5) 



where x is the radius of the sphere. The rms density contrast at scale x is obtained by averaging the above 
expression over all possible locations of the sphere inside the computational box. 



isCg-ikyg-ik'Zg-i(k+k')-ro 



(A6) 



The integral over V^ox reduces to 



(A7) 



We substitute this expression in equation (A6), and use the Kronecker 5 to eliminate the summation over 
k'. Equation (A6) reduces to 



167r2a;<5 



disc j 2 



sph(O) 



The remaining integral can be evaluated easily. Equation (A8) reduces to 

c^l=Y,\6t'?W{kx), 



(A8) 



(A9) 



where 



W{y) = (sin y - y cos yf . 



Using equation (26), we can rewrite this expression in an integral form, 

2 ^bo 



(27r)= 



■ j d^k\5^^''\^W{kx). 



(AlO) 



(All) 
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Table 1. Parameters for the Entire Database 





Ao 


Ho 




n 


^ of dumps 


Codes 




0.20 


0.00 


55 


0.3 


1.2187 


101 


Gbl, 


Gb2, 


Gb3 


0.20 


0.00 


60 


0.3 


1.1539 


92 


Ael 






0.20 


0.00 


65 


0.3 


1.0966 


85 


Acl, 


Ac2, 


Ac3 


0.20 


0.00 


65 


0.5 


1.3188 


85 


Aal, 


Aa2, 


Aa3 


0.20 


0.00 


70 


0.3 


1.0454 


79 


Af 1 






0.20 


0.00 


75 


0.3 


0.9993 


74 


Adl, 


Ad2, 


AdS 


0.20 


0.00 


75 


0.4 


1.1228 


74 


Tal, 


Ta2, 


Ta3 


0.20 


0.00 


75 


0.5 


1.2190 


74 


Abl, 


Ab2, 


Ab3 


0.20 


0.00 


75 


0.6 


1.2979 


74 


Tbl, 


Tb2, 


Tb3 


0.20 


0.00 


75 


0.7 


1.3648 


74 


Tel, 


Tc2, 


Tc3 


0.20 


0.00 


85 


0.3 


0.9191 


65 


Hbl, 


Hb2, 


Hb3 


0.20 


0.80 


55 


0.8 


1.2057 


128 


Gel, 


Gc2, 


Gc3 


0.20 


0.80 


60 


0.6 


0.9948 


117 


Kel 






0.20 


0.80 


65 


0.6 


0.9326 


108 


Kcl, 


Kc2, 


Kc3 


0.20 


0.80 


65 


0.7 


1.0062 


108 


Ual, 


Ua2, 


Ua3 


0.20 


0.80 


65 


0.8 


1.0702 


108 


Kal, 


Ka2, 


Ka3 


0.20 


0.80 


65 


0.9 


1.1269 


108 


Ubl, 


Ub2, 


Ub3 


0.20 


0.80 


65 


1.0 


1.1568 


108 


Ucl, 


Uc2, 


Uc3 


0.20 


0.80 


70 


0.6 


0.8771 


101 


Kf 1 






0.20 


0.80 


75 


0.6 


0.8273 


94 


Kdl, 


Kd2, 


Kd3 


0.20 


0.80 


75 


0.8 


0.9629 


94 


Kbl, 


Kb2, 


Kb3 


0.20 


0.80 


85 


0.8 


0.8749 


83 


Hcl, 


Hc2, 


Hc3 


0.30 


0.00 


75 


0.85 


1.1748 


68 


Eal, 


Ea2, 


Ea3 


0.30 


0.70 


75 


0.9 


0.8796 


81 


Pal, 


Pa2, 


Pa3 


0.35 


0.00 


60 


0.6 


1.0670 


82 


Del 






0.35 


0.00 


65 


0.6 


1.0167 


76 


Dal 






0.35 


0.00 


65 


0.8 


1.1428 


76 


Del 






0.35 


0.00 


70 


0.6 


0.9718 


71 


Df 1 






0.35 


0.00 


75 


0.6 


0.9314 


66 


Dbl 






0.35 


0.00 


75 


0.8 


1.0561 


66 


Ddl 






0.35 


0.65 


60 


0.7 


0.8614 


95 


Mel 






0.35 


0.65 


65 


0.7 


0.8098 


88 


Nal 






0.35 


0.65 


65 


0.9 


0.9251 


88 


Ncl 






0.35 


0.65 


70 


0.7 


0.7640 


82 


Nf 1 






0.35 


0.65 


75 


0.7 


0.7228 


76 


Nbl 






0.35 


0.65 


75 


0.9 


0.8363 


76 


Ndl 






0.50 


0.00 


60 


0.8 


0.9912 


76 


Bel 






0.50 


0.00 


65 


0.8 


0.9457 


70 


Bel, 


Bc2, 


Be3 


0.50 


0.00 


65 


1.0 


1.0439 


70 


Bal, 


Ba2, 


Ba3 


0.50 


0.00 


70 


0.8 


0.9051 


65 


Bf 1 






0.50 


0.00 


75 


0.8 


0.8686 


61 


Bdl, 


Bd2, 


Bd3 
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Table 1 — Continued 





Ao 


Ho 


178 


n 


# of dumps 




Codes 




0.50 


0.00 


75 


1.0 


0.9656 


61 


Bbl, 


Bb2, 


Bb3 


0.50 


0.50 


60 


0.8 


0.8264 


83 


Lei 






0.50 


0.50 


65 


0.8 


0.7808 


77 


Lcl, 


Lc2, 


Lc3 


0.50 


0.50 


65 


1.0 


0.8807 


77 


Lai, 


La2, 


La3 


0.50 


0.50 


70 


0.8 


0.7403 


71 


Lf 1 






0.50 


0.50 


75 


0.8 


0.7049 


66 


Ldl, 


Ld2, 


LdS 


0.50 


0.50 


75 


1.0 


0.8024 


66 


Lbl, 


Lb2, 


Lb3 


0.70 


0.00 


65 


0.9 


0.8461 


64 


Ccl, 


Cc2, 


Cc3 


0.70 


0.00 


65 


1.1 


0.9346 


64 


Cal, 


Ca2, 


Ca3 


0.70 


0.00 


75 


0.9 


0.7773 


56 


Cdl, 


Cd2, 


Cd3 


0.70 


0.00 


75 


1.1 


0.8648 


56 


Cbl, 


Cb2, 


Cb3 


0.70 


0.30 


65 


0.9 


0.7720 


67 


Mcl, 


Mc2, 


Mc3 


0.70 


0.30 


65 


1.1 


0.8601 


67 


Mai, 


Ma2, 


Ma3 


0.70 


0.30 


75 


0.9 


0.7042 


58 


Mdl, 


Md2, 


Md3 


0.70 


0.30 


75 


1.1 


0.7912 


58 


Mbl, 


Mb2, 


Mb3 


1.00 


0.00 


50 


0.5 


0.5836 


75 


Xgl, 


Xg2, 


Xg3 


1.00 


0.00 


55 


1.0 


0.8465 


69 


Gal, 


Ga2, 


Ga3 


1.00 


0.00 


60 


1.0 


0.8057 


63 


Xel 






1.00 


0.00 


65 


0.9 


0.7234 


58 


Sal, 


Sa2, 


Sa3 


1.00 


0.00 


65 


1.0 


0.7698 


58 


Xal, 


Xa2, 


Xa3 


1.00 


0.00 


65 


1.1 


0.8120 


58 


Sbl, 


Sb2, 


Sb3 


1.00 


0.00 


65 


1.2 


0.8506 


58 


Xcl, 


Xc2, 


Xc3 


1.00 


0.00 


65 


1.3 


0.8861 


58 


Scl, 


Sc2, 


Sc3 


1.00 


0.00 


70 


1.0 


0.7380 


54 


Xf 1 






1.00 


0.00 


75 


1.0 


0.7094 


50 


Xbl, 


Xb2, 


Xb3 


1.00 


0.00 


75 


1.2 


0.7893 


50 


Xdl, 


Xd2, 


Xd3 


1.00 


0.00 


85 


1.0 


0.6605 


44 


Hal, 


Ha2, 


Ha3 
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Table 2. Age of the Universe in Gigayears 



On 


An 




1-0 


0.20 


0.00 


55 


15.05 


0.20 


0.00 


60 


13.79 


0.20 


0.00 


65 


12.73 


0.20 


0.00 


70 


11.82 


0.20 


0.00 


75 


11.04 


0.20 


0.00 


85 


9.74 


0.20 


0.80 


55 


19.13 


0.20 


0.80 


60 


17.53 


0.20 


0.80 


65 


16.19 


0.20 


0.80 


70 


15.03 


0.20 


0.80 


75 


14.03 


0.20 


0.80 


85 


12.38 


0.30 


0.00 


75 


10.54 


0.30 


0.70 


75 


12.57 


0.35 


0.00 


60 


12.92 


0.35 


0.00 


65 


11.93 


0.35 


0.00 


70 


11.08 


0.35 


0.00 


75 


10.34 


0.35 


0.65 


60 


15.04 


0.35 


0.65 


65 


13.88 


0.35 


0.65 


70 


12.89 


0.35 


0.65 


75 


12.03 


0.50 


0.00 


60 


12.28 


0.50 


0.00 


65 


11.34 


0.50 


0.00 


70 


10.53 


0.50 


0.00 


75 


9.82 


0.50 


0.50 


60 


13.54 


0.50 


0.50 


65 


12.50 


0.50 


0.50 


70 


11.61 


0.50 


0.50 


75 


10.83 


0.70 


0.00 


65 


10.72 


0.70 


0.00 


75 


9.29 


0.70 


0.30 


65 


11.26 


0.70 


().:-!() 


7.5 


<).7() 


1.00 


0.00 


50 


13.04 


1.00 


0.00 


55 


11.85 


1.00 


0.00 


60 


10.86 


1.00 


0.00 


65 


10.03 


1.00 


0.00 


70 


9.31 


1.00 


0.00 


75 


8.96 


1.00 


0.00 


85 


7.67 
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Table 2 — Continued 
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Table 3. RMS Density Fluctuations from X-ray Clusters 



^0 


Ao 


—clus 


0.20 


0.00 


1.162 


0.30 


0.00 


1.004 


0.35 


0.00 


0.946 


0.50 


0.00 


0.827 


0.70 


0.00 


0.702 


1.00 


0.00 


0.600 


0.20 


0.08 


1.479 


0.30 


0.70 


1.160 


0.35 


0.65 


1.059 


0.50 


0.50 


0.863 


0.70 


0.30 


0.719 


1.00 


0.00 


O.GOO 
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Table 4. Values of for all Simulations 



Oodes 


_num 


Uodes 


—iiuin 


Oodes 




Codes 


_num 
<^8 


Aal 


0.522 


Dal 


0.631 


Kf 1 


0.741 


Sb2 


1.033 


Aa2 


0.575 


Dbl 


0.712 


Lai 


0.903 


Sb3 


1.101 


Aa3 


0.450 


Del 


0.807 


La2 


1.099 


Scl 


1.013 


Abl 


0.579 


Ddl 


0.881 


La3 


1.258 


Sc2 


1.245 


Ab2 


0.479 


Del 


0.606 


Lbl 


0.868 


Sc3 


1.332 


Ab3 


0.475 


Df 1 


0.595 


Lb2 


0.990 


Tal 


0.379 


Acl 


0.269 


Eal 


0.886 


Lb3 


0.966 


Ta2 


0.395 


Ac2 


0.347 


Ea2 


0.902 


Lcl 


0.689 


Ta3 


0.425 


Ac3 


0.280 


Ea3 


0.846 


Lc2 


0.731 


Tbl 


0.600 


Adl 


0.321 


Gal 


1.024 


Lc3 


0.689 


Tb2 


0.618 


Ad2 


0.308 


Ga2 


1.094 


Ldl 


0.843 


Tb3 


0.601 


AdS 


0.289 


Ga3 


0.821 


Ld2 


0.709 


Tel 


0.773 


Ael 


0.268 


Gbl 


0.352 


Ld3 


0.841 


Tc2 


0.694 


Af 1 


0.326 


Gb2 


0.283 


Lei 


0.740 


Tc3 


0.681 


Bal 


0.863 


Gb3 


0.276 


Lf 1 


0.782 


Ual 


0.776 


Ba2 


1.036 


Gel 


0.809 


Mai 


1.145 


Ua2 


0.595 


Ba3 


0.994 


Gc2 


0.749 


Ma2 


1.121 


Ua3 


0.643 


Bbl 


0.828 


Gc3 


0.851 


Ma3 


1.007 


Ubl 


0.707 


Bb2 


1.030 


Hal 


0.940 


Mbl 


1.175 


Ub2 


0.951 


Bb3 


1.056 


Ha2 


0.942 


Mb2 


1.176 


Ub3 


1.003 


Bel 


0.781 


Ha3 


1,053 


Mb3 


1.013 


Ucl 


0.936 


Bc2 


0.889 


Hbl 


0.291 


Mcl 


0.787 


Uc2 


1.195 


Bc3 


0.816 


Hb2 


0.311 


Mc2 


0.797 


Uc3 


0.910 


Bdl 


0.866 


Hb3 


0.281 


Mc3 
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Figure Captions 

Fig. 1. — Left panels: parameter phase space for the database. Models are indicated by dots. The number 
next to each dot indicates the number of simulations. Top panel flo — Xo phase space; middle panel: flo — Hq 
phase space; bottom panel: Oq ~ phase space. Right panels: same as left panels, but for a subregion of 
the parameter phase space. Top panel: Hq = 0.65kms~^Mpc~^ models only; middle and bottom panels: 
Ao = models only 

Fig. 2. — Sample FORTRAN program which reads a file from the database 

Fig. 3. — Relationships between the various values of cts- Each dot represents one simulation. The dashed 
lines indicate the equality between the values of erg plotted. 

Fig. 4. — Two-point correlation function ^ versus separation in Mpc for models with Oq = 0.2 (solid curves). 
The values of the parameters are indicated in each panel. Panels with several curves show ^ for models 
with various values of erg. These values are indicated in the same order as the curves, from top to bottom. 
The dashed line shows the power law given by equation (56). The thick dash indicates the location of the 
softening scale e = 300 kpc. 

Fig. 5. — Same as Fig. 4, except for models with Oq = 0.3 and 0.35. 
Fig. 6. — Same as Fig. 4, except for models with Oq = 0.5 and 0.7. 
Fig. 7. — Same as Fig. 4, except for models with ilo = 1. 
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parameter (n=262144) 
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close(unit=l) 




stop 




end 






"1 — I — I — I — I — I — I — I — I — I — I — I — I — r 



J I I I I I I I L 



J I I L 



.5 



1.5 



disc 



"1 — I — I — I — I — I — I — I — r 



"1 — I — I — r 



. • .1 

I .- : 

I- ! 
1 • 



gauss 



L^-l \ I \ I I I I I I \ I \ L 

.5 1 1.5 

cont 

(^8 



I I I I I I I I I I I I I I I I I I I l_ 
as = 0.3 



o 



Mr 



o 



CD 
O 



o 




no = 0.2 

- Xo = 0.0 
Ho = 70 

I I I I I I I I I 



-I ' 


1 1 1 1 M M M M M M |_ 




CTs = 0.8 


= Oo = 


= 0.2 ^' 


- ^0 = 


= 0.8 ^ 


- H.= 


= 55 I 


1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 h 





M M M M M M M M |_ 




(Tg = 0.6 - 


= Oo = 


= 0.2 I 


~ >^o = 


= 0.8 V 


= Ho = 


= 70 I 


1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 h 


-1 - 


.5 .5 1 1.^ 




log (separation) 



1 1 1,1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ^ 


— 


CTb = 0.3 - 


- 
- 

- % - 


\ 

= 0.2 \. I 


~ Ao = 


= 0.0 \^ 


= Ho = 


= 60 V ; 


n 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 h 



I I I I I I I I I I I I I I I I I I I l_ 
as = 0.5 




_l 1 1 >J 


M M M M M M 1 M M |_ 




CTs = 0.6 


= Oo = 


= 0.2 \-, I 


- Ao = 


= 0.8 ^ 


- Ho = 


= 60 I 


1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 h 



_| 1 'i-j 1 


1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 


'b 






CTs = 0.8 








0.6 




- % - 


= 0.2 






— ^0 = 


= 0.8 






= Ho = 


= 75 






n 1 1 1 1 


1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 


iF 




TTTT 



I I I I I I I I I 



I I I I I I I I I 
CTa = 0.3 - 




no = 0.2 
- Xo = 0.0 
Ho = 85 

I I I I I I I I I 



I I I I 



I I I I 





Mil l 


1 1 1 1 1 1 1 1 1 








CTa = 1.0 








0.9 








0.8 








\ 0.7 








\^ 0.6 




= Oo = 


= 0.2 






- \ = 


= 0.8 






- Ho = 

1 1 1 1 1 


= 65 

Mill 


II 1 II 1 II 1 


iF 



_l i-i, 1 1 


II 1 1 1 


II 1 II 1 II II |_ 






CTs = 0.8 - 


I Oo = 


= 0.2 




~ Ao = 


= 0.8 




= Ho = 


= 85 




1 1 II 1 


1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 h 



-.5 .5 1 
log (separation) 



1.5 -1 -.5 .5 1 
log (separation) 



1.5 



3 
2 
1 



1 

2 
3 

2 

1 



■1 

2 
3 

2 

1 



1 

2 
3 

2 

1 



1 

2 




-.5 .5 1 1.5 
log (separation) 




-1 -.5 .5 1 1.5 -1 -.5 .5 1 1.5 
log (separation) log (separation) 




-1 -.5 .5 1 1.5 -1 -.5 .5 1 1.5 -1 -.5 .5 1 1.5 
log (separation) log (separation) log (separation) 



I I I I I I I I I I I I I I I I I I I l_ 
as = 0.5 




r-i, I I 



I I I I I I I I I I I I I I 

(78= 1-0 




% = 1.0 

- Xo = 0.0 
Ho = 85 

I I I I I I I I I I I I I I I I I I I 



I I I I 



I I I I I I I I I I I I I I I I I I I |_ 

(Ts = 1-0 - 



I I I I I I I I I I I I I I I I I I I |_ 

as = 1.0 = 




no = 1.0 

- Ao = 0.0 
Ho = 70 

I I I I I I I I I I I I I I 



-1 -.5 .5 1 
log (separation) 



1.5 -1 -.5 .5 1 
log (separation) 



-1 -.5 .5 1 
log (separation) 



1.5 



